#
#The Following R-Code creates a set of world-maps which demonstrate the effects of
#our B&P inflation stage covariates, on the probability of non-harmony for our 
#ZiOPC replications of Besley and Persson (2009).

#This corresponds to Figure 1 in the main paper


# clear memory
rm( list=ls() ) 						      

# load some libraries	
library(mvtnorm)
library(car)                                                      
library(scatterplot3d)
library(plotrix)
library(Hmisc)
library(Zelig)
library(foreign, pos=4)
library(pscl)
library(rworldmap)

# set the seed
set.seed(2)

# load dataset
BP.data <- read.dta("bp_exact_for_analysis.dta", 
  convert.dates=TRUE, convert.factors=TRUE, missing.type=TRUE, 
  convert.underscore=TRUE, warn.missing.labels=TRUE)

#summarize dataset
summary(BP.data)

#calculate (and store as variable "NH") the in-sample predicted probabilities of non-harmony 
BP.data$NH<-(pnorm(11.610+(-1.281*BP.data$logGDPpc)+(-0.368*BP.data$parliament)))

#format the predicted probability data for use in the mapping program
dataformap<-with(BP.data, summarize(NH, BP.data[,1], mean, na.rm=TRUE))
dataformap<-as.data.frame(dataformap)
dataformap<-na.omit(dataformap)

#add ISO codes to our predicted probability values for each country included
ISO.2<-rbind("AF", "AL", "DZ","AO", "AR", "AM","AU","AT","AZ","BH","BD","BY","BE","BJ","BO","BA","BW","BR","BG","BF","MM","BI","KH","CM","CA","CV","CF","TD","CL","CN","CO","CR","CI","HR","CU","DK","DJ","DO","EC","EG","SV","GQ","EE","FI","FR","GA","GM","GE","DE","GH","GR","GT","GN","GW","HN","HU","IN","ID","IR","IQ","IE","IL","IT","JM","JP","JO","KZ","KE","KW","KG","LA","LV","LB","LS","LR","LY","LT","MK","MG","MW","MY","ML","MR","MU","MX","MD","MN","MA","MZ","NA","NP","NL","NZ","NI","NE","NG","KP","NO","OM","PK","PA","PY","PE","PH","PL","PT","QA","RO","RU","RW","SA","SN","RS","SC","SL","SG","SK","SI","SO","ZA","KR","ES","LK","SD","SZ","SE","CH","SY","TJ","TZ","TH","TG","TT","TN","TR","TM","UG","UA","AE","GB","US","UY","UZ","VE","VN","YE","Y3","CD","ZM","ZW") 
dataformap<-cbind(dataformap,ISO.2)
dataformap<-as.data.frame(dataformap)

#generate values for mapping
sPDF<-joinCountryData2Map(dataformap
     , joinCode = "ISO2" 
     , nameJoinColumn = "ISO.2"
     , projection="none"  
     )

#Create a colored map for Pr(Non-Harmony)
mapParams <- mapCountryData( sPDF, nameColumnToPlot="NH", addLegend=FALSE, mapTitle="")
do.call(addMapLegend, c(mapParams,legendLabels="limits",legendWidth=0.5,legendIntervals="data"))

#Create a black and white map for Pr(Non-Harmony)
mapParams <- mapCountryData( sPDF, nameColumnToPlot="NH", colourPalette="white2Black", addLegend=FALSE, mapTitle="")
do.call(addMapLegend, c(mapParams,legendLabels="limits",legendWidth=0.5,legendIntervals="data"))

